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Abstract 

This  paper  examines  uniqueness  and  stability  results  for  an  inverse  problem  in  thermal 
imaging.  The  goal  is  to  identify  an  unknown  boundary  of  an  object  by  applying  a  heat  flux 
and  measuring  the  induced  temperature  on  the  boundary  of  the  sample.  The  problem  is 
studied  both  in  the  case  in  which  one  has  data  at  every  point  on  the  boundary  of  the  region 
and  the  case  in  which  only  finitely  many  measurements  are  available.  An  inversion  procedure 
is  developed  and  used  to  study  the  stability  of  the  inverse  problem  for  various  experimental 
configurations. 
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1  Introduction 


Thermal  imaging  is  a  technique  of  wide  utility  in  non-destructive  testing  and  evaluation.  The 
technique  is  used  to  recover  information  about  the  internal  condition  of  an  object  by  applying 
a  heat  flux  to  its  boundary  and  observing  the  resulting  temperature  response  on  the  object’s 
surface.  From  this  information  one  attempts  to  determine  the  internal  thermal  properties  of 
the  object,  of  the  shape  of  some  unknown  portion  of  the  ])oundary.  'Thermal  imaging  has 
been  much  investigated  as  a  method  for  detecting  damage  or  corrosion  in  aircraft.  See  [10] 
for  an  account  of  the  technology  and  typical  data  processing  techniques  that  are  employed, 
and  a  more  extensive  bibliography  on  the  subject. 

One  of  the  most  common  uses  of  thermal  imaging  is  for  the  detection  so-called  “back 
surface”  corrosion  and  damage.  Briefly,  one  attempts  to  determine  whether  some  inaccessible 
portion  of  an  object’s  boundary  has  corroded,  and  therefore  changed  shape.  In  this  paper 
we  investigate  the  inverse  problem  of  determining  changes  in  the  boundary  proflle  of  a  two- 
dimensional  sample  by  using  thermal  imaging.  We  consider  a  certain  portion  of  the  surface 
of  a  rectangular  sample  to  be  accessible  for  measurements  and  the  remainder  of  the  surface, 
which  may  be  corroded,  inaccessible.  This  problem  has  been  considered  by  others  [3,  4]  with 
an  emphasis  on  recovering  estimates  of  the  unknown  surface  from  data  by  using  an  output 
least-squares  method. 

We  examine  both  a  continuous  and  finite  data  version  of  the  inverse  problem.  The 
continuous  version  assumes  that  one  has  data  at  every  point  on  the  accessible  portion  of  the 
object’s  surface.  The  finite  data  version  assumes  that  only  finitely  many  measurements  have 
been  made.  Our  goals  are 

•  To  examine  uniqueness  and  continuous  dependence  results  for  the  continuous  version 
of  the  inverse  problem,  and  what  they  imply  for  the  finite  data  inverse  problem. 

•  To  examine  how  various  experimental  parameters  affect  stability  and  resolution  for  the 
finite  data  inverse  problem,  especially  the  eflfect  of  measurement  locations  on  stability. 

•  To  determine  how  one  might  incorporate  a  priori  information  or  assumptions  into  the 
finite  data  inverse  problem. 
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Our  main  focus  is  not  to  develop  inversion  algorithms,  but  in  the  course  of  examining  the 
problem  we  derive  (but  do  not  prove  convergence  for)  an  inversion  procedure  for  the  finite 
data  inverse  problem.  This  algorithm  allows  the  easy  incorporation  of  a  priori  assumptions 
into  the  inversion  process.  We  apply  the  algorithm  to  several  simulated  data  sets  to  illustrate 
our  conclusions.  Our  study  of  the  stability  of  the  inverse  problem  reduces  to  studying  the 
invertibility  of  a  certain  matrix,  which  we  do  with  a  singular  value  decomposition.  We  do 
not  make  any  explicit  finite  dimensional  parameterization  of  the  unknown  surface. 

We  should  note  that  a  very  similar  approach  has  been  used  in  [8]  to  study  resolution  and 
stability  for  the  inverse  conductivity  problem.  Isaacson,  Cheney,  and  others  [6,  7]  have  also 
carried  out  similar  sensitivity  studies  related  to  the  inverse  conductivity  problem,  especially 
the  effect  of  finitely  many  measurements  on  the  inversion  process. 

The  outline  of  the  paper  is  as  follows.  In  Section  2  we  present  the  mathematical  for¬ 
mulation  of  the  continuous  and  finite  data  versions  of  the  inverse  problem.  In  Section  3  we 
derive  a  linearized  version  of  the  continuous  problem,  and  in  Section  4  we  show  how  this 
leads  (as  thermal  inverse  problems  often  do)  to  a  first  kind  integral  equation  which  must  be 
inverted.  In  Section  5  we  use  the  integral  equation  formulation  to  examine  uniqueness  and 
stability  results  for  the  linearized  version  of  the  inverse  problem.  In  Section  6  we  consider  an 
algorithm  for  solving  the  finite  data  version  of  the  inverse  problem  and  how  this  approach 
can  be  used  quantify  the  stability  of  the  problem.  In  the  last  section  we  present  a  variety 
of  numerical  studies  to  examine  the  effects  that  various  experimental  parameters  have  on 
the  stability  and  resolution  of  the  inversion  process,  and  the  effect  of  incorporating  a  priori 
assumptions  into  the  inversion  procedure. 


2  The  Inverse  Problem 

Consider  a  sample  to  be  imaged  as  a  two-dimensional  region  lying  between  the  two  surfaces 
;r2  =  S{xi)  and  ;r2  =  1  as  illustrated  below. 
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_  X2=S(Xi) 

X 

1 

Figure  1:  Sample  geometry. 

We  will  refer  to  the  surface  X2  =  1  as  the  “top”  or  “front”  surface  and  X2  =  S{xi}  as  the 
“back”  surface.  We  assume  that  the  ends  of  the  sample  are  sufficiently  far  away  that  they 
can  be  ignored,  so  for  our  purposes  the  sample  is  unbounded  in  the  .ti  direction.  The  top 
surface  is  assumed  to  be  accessible  for  inspection  or  measurements,  but  the  back  surface 
X2  =  S{xi)  is  inaccessible.  This  is  the  portion  of  the  sample  to  be  inspected  for  corrosion. 
The  ideal  uncorroded  case  is  a  flat  back  surface  S{xi)  =  0.  In  the  corroded  case  above 
5(.ti)  >  0  for  some  values  of  .ti.  We  will  assume  that  the  function  S  belongs  to  H^(JR), 
although  this  assumption  will  later  be  relaxed.  In  particular,  since  H^{]R)  C  C^(]R)  there  is 
a  continuous  unit  normal  vector  field  on  the  back  surface.  The  goal  is  to  determine  the  back 
surface  or  the  function  S  by  taking  measurements  only  on  the  front  surface. 

A  time-dependent  heat  flux  is  applied  to  the  top  of  the  sample  X2  =  1.  We  will 

assume  that  the  sample  material  is  homogeneous  with  thermal  difirusivity  k  and  thermal  con¬ 
ductivity  o,  both  known  constants.  We  will  use  T{x,t)  to  denote  the  resulting  temperature 
induced  in  f2,  where  ;r  =  (;ri,.T2).  The  direct  thermal  diffusion  problem  will  be  modeled  by 
the  standard  heat  equation 

=  0  in  O,  (2.1) 

=  on  .T2  =  1, 

=  0  on  X2  =  S(xi), 

=  To(x), 

for  #  >  0,  where  ^  denotes  the  outward  normal  derivative  on  the  boundary  of  Q.  The 
function  To(.r)  is  the  initial  temperature  of  the  region  ft  at  time  t  =  0.  Note  that  the  back 
surface  is  assumed  to  block  all  heat  conduction. 


dt 


-K  AT 


dT 


a- 


du 

du 
T(.t,0) 


a- 
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We  will  also  assume  that  the  heat  flux  g{xi,t)  is  periodic,  of  the  form  Re[.g(.Ti)e*‘^'*]  with 
LO  >  0.  For  simplicity,  we  also  take  the  constants  k  and  a  equal  to  one.  Under  these 
assumptions  the  solution  to  equation  (2.1)  is  given  as  T{x,t)  =  Re[e’^’U<(.r)]  where  u{x) 
satisfies 


An  —  iixu 
dll 

dll 


0  in  Q, 

g{xi)  on  X2  =  1, 
0  on  Xo  =  5(:ri), 


(2.2) 


at  least  after  transients  from  the  initial  condition  have  died  out.  The  main  case  of  interest  is 
that  in  which  g{x\)  is  constant,  corresponding  to  uniform  heating  of  the  outer  surface.  This 
is  typically  the  case  when  heat  or  flash  lamps  are  used  to  provide  the  input  flux  g.  For  the 
moment,  however,  we  will  not  restrict  g. 

There  are  two  versions  of  the  inver.se  problem  to  be  considered: 


Continuous  Version: 

Given  measurements  of  n(.T)  at  all  points  on  the  top  surface  .T2  =  1,  determine  the 
function  S{xi). 

Finite  Data  Version: 

Given  measurements  of  ii{x)  on  the  top  surface  X2  =  1  at  points  .Ti  =  Oi,G2,.  . .  ,an, 
estimate  the  function  5(;ri). 

The  finite  data  version  corresponds  to  the  case  in  which  one  has  actual  measurements.  The 
data  need  not  be  actual  point  measurements  of  the  temperature  u,  but  this  is  the  most 
common  situation.  Of  particular  interest  are  the  questions 

1.  Can  the  function  S{xi)  be  uniquely  determined  by  knowing  u{xi)  for  all  Xi  on  the  top 
surface? 

2.  If  5(;ri)  is  uniquely  determined  by  w(.Ti),  how  sensitive  is  S{xi)  to  perturbations  in  the 
data?  What  kinds  of  features  in  the  back  surface  .T2  =  S'(.ti)  can  or  cannot  be  easily 
determined  from  the  data? 
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3.  Since  any  practical  application  falls  under  the  finite  data  formulation,  how  stable  is 
the  estimate  of  S{xi)  based  on  finitely  many  pieces  of  data?  What  factors  influence 
stability  in  this  case,  and  is  there  an  inversion  procedure  to  produce  a  reasonable 
estimate  of  5(;ri)  using  finitely  many  measurements? 

The  first  question  is  easily  answered  “yes”  by  a  standard  argument.  This  is  the  content  of 
the  following -result. 

Theorem  2.1  (Uniqueness)  Let  u(xi,X2‘,S)  denote  the  solution  to  (2.2)  with  back  surface 
S  and  nonzero  flux  g.  If  u{xi.,l\Si)  =  u(;ri,l;52)  for  each  (;ri,l)  in  an  open  subset  C  of 
the  top  surface  of  Cl,  then  Si  =  82- 

Proof:  Suppose  ^  82-  Using  the  shorthand  notation  u,-  =  u{xi,X2',  Si),  we  have  that  ui 
and  U2  have  the  same  Cauchy  data  on  C,  and  by  unique  continuation  Ui  and  U2  agree  on  any 
connected  domain  on  which  both  are  defined,  provided  that  that  domain  contains  an  open 
portion  of  C.  Assume  that  Si  and  S2  are  not  equal  and  denote  by  Clj  the  region  bounded 
by  X2  =  1  and  X2  =  Si{xi).  Si  7^  S2,  so  there  is  some  non-empty  connected  component 
D  of  Qi\  CI2  or  ^2  \  Let  us  suppose  the  latter  so  that  the  region  D  is  bounded  by 
X2  =  5i(xi)  and  X2  =  S2{xi).  On  X2  =  52{.ti)  we  know  that  the  normal  derivative  of  U2  is 
identically  zero;  on  X2  =  S'i(;ri)  we  know  that  the  normal  derivative  of  tq  is  zero  (from  inside 
rii)  and  since  U2  =  Ui  and  U2  is  smooth  across  X2  =  5'i(.Ti),  we  conclude  that  the  normal 
derivative  of  U2  vanishes  on  the  boundary  of  D.  This  forces  U2  =  0  inside  D.  Standard  elliptic 
regularity  argiiments  force  112  =  0  inside  ^27  which  implies  that  the  flux  g  is  identically  zero, 
a  contradiction.  Thus  we  must  have  =  S2  and  so  the  back  surface  5(xi)  is  uniquely 
determined  by  the  boundary  data  on  any  open  portion  of  the  top  surface.  ■ 

The  second  and  third  questions  will  be  examined  in  the  next  few  sections  by  considering 
a  linearization  of  the  original  problem. 


3  A  Linearization 

We  now^  linearize  the  original  direct  problem  given  by  equation  (2.2)  with  respect  to  the 
function  S.  Let  uq{x)  denote  the  solution  to  (2.2)  with  5  =  0.  The  surface  X2  =  0  is  the 
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point  at  which  we  will  linearize,  since  this  represents  the  imcorroded  or  ideal  profile  from 
which  we  hope  to  detect  any  deviation. 

Let  flo  denote  the  region  {(;ri,.r2)  :  0  <  .T2  <  1}  in  For  any  S  G  //'^(R)  C  C^(R), 
let  us  construct  the  map  0  from  ft  to  Hq  by 

=  ). 

If  |5(.ri)|  <  1  then  it  is  easy  to  check  that  d>  is  invertible  on  flo-  This  map  fixes  the  top 
boundary  of  and  maps  the  bottom  surface  to  .r2  =  0.  Let  y  =  0(.r)  and  v(y)  =  u{(l)~^{y))  = 
u(x).  Under  such  a  change  of  coordinates  Vj.  =  {D(f))^Vy  and  ^  ^  where  7]y  =  {D(t))u^, 

{D(j))  is  the  derivative  of  d,  and  is  the  transpose.  Under  this  change  of  coordinates 

the  boundary  value  problem  (2.2)  becomes 

V  •  kVv  —  iujv  =  0  in  Uq,  (3.3) 

—  =  g{xi)  on  .r2  =  1, 
dv 

—  =  0  on  X2  =  0, 

07] 

where  k  =  {D4>){D(}))'^ .  The  vector  t]  can  be  written  as 


on  .r2  =  1 .  7]  = 


l-S(xi) 


1  5'(xi) 

''1  +  (S'(.Ti))2  1  +  (S'(xi))^ 

5(xi)-l 


on  .T2  =  0. 


This  change  of  coordinates  removes  the  unknown  S  from  the  definition  of  the  boundary  of 
Q  and  puts  S  into  the  coefficients  of  the  heat  operator  and  boundary  conditions. 

Now  we  linearize  the  problem  with  respect  to  S  about  5  =  0  by  assuming  that  S  =  eS 
for  some  function  S,  where  e  is  some  small  real  number.  Let  w  denote  the  solution  to  (3.3) 
for  a  general  S  and  Uq  the  solution  to  equation  (3.3)  for  5  =  0.  Suppose  that  the  resulting 
perturbation  in  t<o  can  be  represented  as  w  =  Uq  +  ci’o-  If  we  substitute  these  relations  into 
(3.3),  use  the  fact  that  vq  satisfies  (3.3)  with  5  =  0,  and  drop  all  quadratic  e  terms  then  we 
obtain  a  linearized  version  of  the  direct  problem, 

Ai'o  —  iiovQ  =  —  V  ■  [kVwo],  in  flo,  (3-4) 

^  =  -S{yi)g{yi),  on  ^2  =  1 
^  =  -5'(?/i)wo(yi),  on  2/2  =  0, 
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where  k{y)  is  the  matrix 


0  {y2  -  ^)S'{yi) 

k{y)  = 

_(t/2-imyi)  25{yi) 

In  particular,  w'e  are  interested  in  the  special  case  g  =  \  corresponding  to  spatially 
uniform  periodic  heating  of  the  top  surface.  Under  this  condition  the  function  uq  depends 
only  on  X2  and  equation  (3.3)  for  uq  becomes  a  two-point  boundary  value  problem 

Uq  —  iojUQ  =  0  on  (0,1),  (3.5) 

«i(0)  =  0, 

«;{!)  =  1. 

If  we  now  use  u  rather  than  vq  to  denote  the  temperature  perturbation  satisfying  equation 
(3.3)  then  the  linearized  version  of  the  problem  for  g  =1  is  then 

Au  —  iwu  =  {1  —  X2)uq(x2)S''{xi) —  2uq{xi)S{xi)  on  fio?  (3-6) 

^  =  -S{xi)  on  X2  =  1, 
av 

du  ^  ^ 

—  =  0  on  X2  =  0. 
av 

For  simplicity,  this  is  the  version  of  the  problem  we  will  examine,  although  the  more 
general  linearized  version  (3.4)  can  be  examined  using  similar  techniques.  Note  that  the 
linearized  problem  is  defined  on  the  domain  Oq  which  does  not  depend  on  5. 

4  An  Integral  Identity 

Let  the  function  d{a)  =  u{a,  1)  denote  the  top  surface  data  from  the  direct  problem  (3.6). 
Given  that  the  relation  between  S  and  d  is  linear,  it  seems  reasonable  that  this  relationship 
can  be  expressed  by  an  integral  operator 

/OO 

My)S{y)  dy 

“OO 

where  da  is  some  function  which  depends  on  the  particular  point  a.  Such  a  relationship  does 
exist  and  we  can  say  quite  a  bit  about  the  function(s)  da-,  as  we  now  demonstrate. 
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Let  L  =  A  —  iuj.  By  Green’s  second  identity 

f  f  (  dll  dd\  ,  f  ( ,du  d4)\ 

-  iiL0)  d.v  =  ^  [6-  -  u-j  d.r,  +  (^0-  -  n-j  d.y 

where  0  is  any  sufficiently  regular  function  defined  on  fl.  Assume  that  u  is  a  solution  to 
(3.6)  and  that  0  is  a  function  which  satisfies  L0  =  0  on  flo  and  |^  =  0  on  the  surface  X2  =  0. 
Then  the  above  equation  becomes 

/  c^(;r)((l-.T2)»'o(:r2)5''(:ri)-2u"(.ri)5'(.ri))d.T  +  T  0(.ri,  l)5(.ri)  dxi 

■/f!o  J-oo 

Let  us  now  complete  the  specification  of  (f)  by  requiring  on  the  top  surface,  where  5a 

denotes  a  delta  function  on  the  surface  X2  =  1  at  the  point  .Ti  =  a.  If  we  write  da  to  denote 
the  dependence  of  d  on  a  then  we  obtain 

/  da{-v){{l-X2)u'o{x2)S"{xi)-2ii'^{xi)S{xi))dx+  [  daixiA)S{xi)dxi  = -d{a).  (4.7) 
Jflo  J-oo 

Note  that  this  equation  involves  no  unknown  quantities  except  5  on  the  left  side. 

It  is  worth  saying  a  few  words  about  the  function  da  which  satisfies 


Ada  —  ix)da  =  0  in  Hq, 

dda  j.  T 

—  =  5a  on  ;r2  =  1, 

OP 

dda  f. 

^r—  =  0  on  .V9 

OP 


=  0  on  Xo  =  0. 


Let  r(.T)  be  a  Green’s  function  for  the  operator  A  —  iuj  on  IR^;  such  a  function  is  given  by 

r(.r)  =  — ^(ker(rv/u’)  +  tkei(r\/u’)) , 

aTT 

where  r  =  |.r|  =  ^Jxl  +  .rl  and  ker()  and  kei()  are  the  Kelvin  functions  (see  [1],  section  9.9). 
The  function  ker(7’)  has  a  —  ln(?’)  singularity  as  r  approaches  zero,  while  kei(7')  is  bounded. 
Both  functions  and  their  derivatives  are  smooth  away  from  zero  and  rapidly  decreasing  as  r 
tends  to  infinity,  where  “rapidly  decreasing”  means  faster  than  any  power  of  If  we  define 
ipa{x)  =  — 2r(;r  —  .Tq)  where  Xa  is  the  point  (a,  1)  on  the  top  surface,  then  standard  potential 
theory  arguments  ([9],  chapter  3)  show  that 

Ada  —  i^da  = 
dda  _ 
dp 
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0  in  fio, 

5a  on  X2  =  1. 


(4.9) 


It  is  not  true  that 


=  0  on  X2  =  0. 


However,  if  we  take  Va  G  to  satisfy 


Av„  —  iu)v„ 


0  in  Oo, 


(4.10) 


—  =  0  on  .T2  =  1, 
oiy 

dVa  di’a  „ 

—  =  — —  on  X2  =  0. 

ou  au 

then  (j)a  =  i'a  +  t’a-  Since  the  Neumann  data  for  Va  on  the  bottom  surface  is  —  ^  which 
is  in  (the  singularity  for  ipa  bes  on  the  top  surface,  and  away  from  this  singularity 

ipa  is  smooth  and  rapidly  decreasing),  one  can  show  that  Va  is  in  As  a  result,  the 

function  (pa{x)  has  a  — Mn  |.t|  singularity  near  .r  =  0  and  otherwise  is  smooth  and  rapidly 
decreasing  in  |;r|,  along  with  its  derivatives  of  all  orders. 

If  we  write  the  integrals  in  equation  (4.7)  with  limits,  we  find  that  S  must  satisfy 

/OO 

{Pa{xi)S"{xi)  +  qa{xi)S{xi))  dxi  =  -d{a) 

-OO 

where 

Pa(.ri)  =  /  0a(-Ti,*T2)«o(-'^'2)(l -3:2)d;r2,  (4.11) 

Jo 

=  -2  /  <^a(-Tl,.T2)Mo(.T2)d^:2  +  0a(^l,l)-  (4.12) 

J  0 


(4.11) 

(4.12) 


One  can  check  that  the  integral  in  q^xi)  is  continuous  as  a  function  of  Xi,  smooth  away 
from  .Ti  =  a,  and  rapidly  decreasing  in  xi.  Also,  since  (pa{xi,  1)  has  a  logarithmic  singularity 
at  ;ri  =  a,  so  does  qa-  Moreover  qa  is  an  function.  The  function  Pai-'Pi)  is  also  clearly 
smooth  away  from  xi  =  a  and  rapidly  decreasing  in  xi.  The  singularity  of  looks 

like  the  singularity  of  ker|.T  —  Xa|,  and  one  can  use  this  fact  to  expand  the  function  d'a(.Ti,  .T2) 
near  the  singularity  to  show  that  the  function  Pai^i)  is  actually  in  //^(IR).  Since  both  Pa{xi) 
and  9a(a’i)  tend  rapidly  to  zero  as  |xi[  ->  00  we  can  integrate  by  parts  twice  to  find  that 


/OO  TOO 

Pa (.r  1 ) S" (.Ti )dxi=  /  pI {xi)S (.r  1 )  d.-r  1 . 

■OO  J  —  00 


Equation  (4.9)  can  now  be  written 


/OO 

c„(.ti)5(.ti)  dxi  =  d{a) 

-CO 


(4.13) 
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where 

c.(.ti)  =  -p'Hi-i)  -  €  I"(IR).  (4.14) 

Given  the  translation  invariance  of  this  problem  in  the  .Ti  direction  and  the  fact  that  the 
flux  g{.vi)  =  1  does  not  depend  on  .Ti,  it  is  clear  that  Ca(.ri)  =  c(.ti  —  o)  for  some  function  c. 

Lemma  4.1  The  function  c(.ri)  =  iio{0)ruj(t>o(xi,0) ,  where  do  satisfies  (4-8)  with  a  =  0  and 
uo  satisfies  (3.5). 

In  proving  this  result,  We  will  make  use  of  the  following  simple  fact. 

Remark  4.1  If  f{t)  and  g{t)  are  functions  defined  on  [0,1]  with  /'(O)  =  f'{l)  =  0  and 
^'(0)  =0,  ^'(1)  =  1,  and  g”  =  iiog  then 

C  f"{t)g'{t){l  -t)dt  =  /(I)  -  2ia;  ['  f{t)g(t)  dt  +  iujf{0)g{0)  +  iu  l\l  -  t)f(t)g'(t)  dt. 
.10  Jo  Jo 

To  see  this  identity  requires  only  a  few  applications  of  integration  by  parts,  /  udv  =  uv  — 
f  vdu.  Take  u  =  g'{t){l  —  t),  dv  =  /"(/)  dt  to  obtain 

r  f{t)g'it){l  -t)dt=  ['  f{t)g'{t)  dt  -  iuJ  f  f'{t)g{t){l  -  t)  dt, 

Jo  Jo  Jo 

where  we  have  made  use  of  g"  =  iug.  Integrate  the  first  integral  by  parts  with  u  =  g'  and 
dv  =  f  dt:  integrate  the  second  integral  by  parts  with  u  —  {1  —  t)g  and  dv  =  f  dt.  After  a 
number  of  cancellations,  the  Remark  follows  immediately. 

Proof  of  Lemma  4.1:  From  equations  (4.11),  (4.12)  and  (4.14)  we  obtain 

fl  52  a 

(’((^i)  =  -  ((1  -  .T2)^p^Wo(.r2)  -  2doUo{xi))  d.V2  -  doi^i^  !)• 

J  0  C/X  j 

Recall  that  t/,o  is  a  function  of  .T2  only,  and  that  Uq  =  iuuo  with  Uo(0)  =  0  and  Uo(l)  =  1. 
Since  Ado  —  io^do  =  0  we  have  =  ieodo  —  With  this  substitution  and  Uq  =  ivouo 

c{xi)  =  -  [?w(l  -  X2)doUQ{x2)  -  (1  -  x2)uo^-t  “  dx^  -  0o(2q,  1). 

Jo  0X2 

To  finish  the  proof,  use  Remark  4.1  above  with  f  =  do  and  g  —  uq-  After  a  few  cancellations 
the  statement  of  the  lemma  follows.B 
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Lemma  4.1  implies  a  number  of  things  about  the  function  c(.ri),  in  particular,  that  c  is 
in  It  is  important  to  note  that  t<o(0)  is  never  zero,  for  uq  satisfies  the  two-point 

boundary  value  problem  (3.5)  and  one  easily  finds  that 


uo{x2)  = 


eCCX2  ^  g-a.T2 

a:(e"  -  e-“)  ’ 


where  q:  =  (1-1-  In  particular. 


*^o(0)  =  - r. 

a:(e“  -  e-“) 

The  numerator  above  is  never  zero,  hence  tfo(O)  7^  0.  As  a  result  the  function  c(.Ti)  will  not 
be  identically  zero  for  any  a;  >  0. 

There  is  one  more  fact  which  will  be  extremely  useful.  Defining  the  Fourier  transform 
f{y)  of  a  function  f{x)  of  a  single  variable  by 


A  /*oo 

f{y)  =  /  dx, 

J—oo 


we  have 


Lemma  4.2  For  c  given  by  Lemma  4-i, 

_  ,  2^a;^^o(0) 

c  y)  =  -r-r - 

a(e"  —  e 

where  a  =  y/y^  +  loj.  Also,  the  function  c{y)  is  never  equal  to  zero. 


where  a  =  y/y^  +  ioj.  Also,  the  function  c{y)  is  never  equal  to  zero. 

Proof:  Let  (i>o{yi,X2)  denote  the  Fourier  transform  of  (l>o{xi,X2)  with  respect  to  Xi  only, 
where  do  satisfies  equation  (4.8)  with  a  =  0.  Then 


<t>o{yi,X2)=  /  (i>o{xi,X2)e~^''^y^  dxi- 

J—oo 


If  do  is  sufficiently  rapidly  decreasing  (as  in  our  case)  then  ^  and  the  Fourier  transform 
operator  commute.  Fourier  transforming  both  sides  of  the  boundary  value  problem  (4.8)  with 
respect  to  .ti  yields  a  two-point  boundary  value  problem  for  do{yi,X2)  in  the  X2  variable, 

-  {yl  +  iuj)do  =  0, 

>  = 


■(0)  =  0, 


(4.15) 
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where  we  have  used  5o  =  1.  The  solution  is 


0o(?/ii^'2)  — 


^aT2  ^  g,— 0X2 


o(e"  —  e“°) 

where  a  =  +  iuj.  On  the  bottom  surface  X2  =  0, 

2 

^o(yi,0)  =  — — - 

o(e^  -  e  “) 

which,  combined  with  Lemma  4.1  yields  the  conclusion.  The  fact  that  c{y)  is  never  zero 
follows  from  w  >  0  and  ko(0)  7^  0.  ■ 


5  The  Linearized  Inverse  Problem 

In  this  section  we  will  examine  the  continuous  version  of  the  linearized  inverse  problem. 
Suppose  that  we  have  top  surface  data  d{xi)  —  «(;ri,l)  where  ii  satisfies  equation  (3.6). 
Based  on  equation  (4.13)  and  the  above  noted  fact  that  Ca(;ri)  =  c(.Ti  —  a)  we  conclude  that 
S  and  d  satisfy  the  relationship 

roo 

/  S{yi)c{yi-xi)dyi  =  d{xi) 

J  —  00 

or 

[  S(;f/i)d>(;ri  -  yi)dyi-S*6  =  d{xi)  (5.16) 

J  — 00 

where  (h{xi)  =  c(— ;ri)  and  denotes  convolution.  With  d(xi)  considered  known  this 
becomes  a  first  kind  integral  equation  for  the  unknown  function  5.  The  kernel  (f)  is  known, 
or  can  be  determined,  according  to  Lemma  4.1.  First  kind  integral  equations  have  been 
extensively  studied  ([11,  12])  and  are  well-known  to  be  unstable;  small  perturbations  in  the 
right  hand  side  d{x)  can  lead  to  arbitrarily  large  changes  in  the  solution  5.  However  this 
formulation  of  the  inverse  problem  as  an  integral  equation  will  allow  us  to  obtain  uniqueness 
and  stability  estimates  for  the  linearized  version  of  the  inverse  problem. 

Theorem  5.1  (Uniqueness)  If  the  data  d(x)  is  an  function  and  if  there  exists  a  solution 
S  6  L^(IR)  to  the  linearized  inverse  problem,  then  S  is  unique. 
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Proof:  Suppose  Si  and  S2  are  functions  which  both  give  rise  to  data  d.  Let  S  =  Si  —  52. 
Linearity  implies  that  the  data  for  5  is  identically  zero.  Fourier  transform  both  sides  of  the 
equation  S  *  (j)  =  0.  By  the  convolution  theorem  and  the  fact  that  0  =  c  if  (p{x)  =  c(— .t) 
(where  the  overbar  denotes  complex  conjugation)  we  obtain  ^S  =  0.  By  Lemma  4.2  d  = 
c  7^  0  and  we  conclude  that 

5  =  0. 

so  that  5  =  0  or  5i  =  52.  ■ 

Remark  5.1  In  the  preceding  proof,  we  assumed  a  priori  that  5  is  in  L^(1R.).  In  general, 
for  an  arbitrary  d  G  L^(IR)  we  cannot  find. a  function  S  in  which  gives  rise  to  data  d  via 
equation  (5.16). 

The  convolution  equation  (5.16)  also  provides  information  on  continuous  dependence.  Since 
the  function  is  smooth  and  never  equal  to  zero,  we  can  define  the  space  of  functions  L*(]R) 
with  the  norm 


From  Lemma  4.2  it  follows  that  grows  like  ze^.  The  norm  ||  ||*  thus  puts  a  heavy 
penalty  on  high  frequencies;  the  functions  in  this  space  are  very  smooth.  We  can  Fourier 
transform  both  sides  of  equation  (5.16),  divide  by  ^0  and  take  the  norms  of  both  sides  to 
obtain 

Theorem  5.2  (Continuous  Dependence)  If  a  back  surface  X2  =  S(xi)  generates  front  sur¬ 
face  data  d(x)  for  the  linearized  problem  (3.6),  then 

l|S||i>  <  C||rf||. 

where  C  is  independent  of  d. 

Estimates  of  5  from  data  d  will  thus  be  extremely  sensitive  to  any  noise,  because  the  in¬ 
version  process  weights  a  frequency  /  in  the  data  by  a  factor  proportional  to  fe^.  Lemma  4.2 
and  the  structure  of  the  convolution  operator  mapping  5  to  the  data  d  make  it  clear  that  it 
will  be  difficult  to  estimate  the  high  spatial  frequency  components  in  the  Fourier  decompo¬ 
sition  of  5,  for  these  components  are  heavily  damped  out  by  the  forward  mapping. 
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6  The  Case  of  Finitely  Many  Measurements 


Suppose  that  we  have  point  estimates  rf(a,)  =  ?/(«,  .  1)  of  the  temperature  on  the  top  surface 
at  n  distinct  points.  How  can  we  construct  a  reasonable  estimate  of  the  function  5(;ci)? 
How  can  we  quantify  the  stability  of  the  reconstruction  with  respect  to  errors  in  the  data, 
and  how  does  the  choice  of  measurement  locations  Oj-  affect  the  stability?  Let  us  assume 
that  we  seek  an  estimate  S  G  L^(IR).  Physical  considerations  make  it  desirable  to  obtain  an 
estimate  with  more  regularity,  but  this  will  be  a  consequence  of  the  proposed  reconstruction 
procedure.  Based  on  the  convolution  equation  (5.16)  we  know  that  S  must  satisfy  the  n 
constraints 

/OO 

S{xi)c,{xi)dxi  =  d{ai),  7'  =  l,...,n,  (6.17) 

-OO 

with  c,(.Ti)  =  c(a,  —  Xi)  where  c(.ri)  is  the  function  from  Lemma  4.1  and  <  f,g  >=  fji  fg  is 
the  usual  L‘^  inner  product.  Note  that  since  c,  is  an  function,  S  i-4<  S',  c,  >  is  a  bounded 
linear  functional  on  The  set  (6.17)  is  a  horribly  underdetermined  set  of  equations.  We 
can  expect  to  find  an  entire  translated  subspace  of  functions  of  codimension  n  in  L‘^(1R) 
which  satisfy  the  given  conditions,  and  any  such  function  “solves”  the  inverse  problem,  in 
the  sense  that  it  gives  rise  to  the  measured  data. 

One  practical  method  for  specifying  a  unique  function  in  which  solves  the  inverse 
problem  is  to  seek  that  element  in  which  satisfies  the  given  conditions  and  has  minimal 
norm.  That  such  an  element  exists  follows  from  the  fact  that  the  relations  (6.17)  define  a 
closed  convex  subset  of  and  hence  this  subset  has  a  unique  element  of  minimal  norm.  This 
idea  has  been  used  before  ([8])  to  construct  a  “pseudo-inverse”  for  the  finite  measurement 
case  and  to  characterize  the  stability  and  information  content  for  the  inverse  conductivity 
problem,  and  has  also  been  used  for  reconstruction  from  partial  information  in  tomographic 
problems  [5].  The  approach  has  several  merits:  In  the  present  case  it  leads  to  an  exceptionally 
easy  and  efficient  inversion  algorithm  which  allows  us  to  stud}’  the  conditioning  of  the  inverse 
problem  independent  of  any  explicit  finite  dimensional  parameterization  of  the  unknown  S. 
By  weighting  the  space  appropriately  we  can  also  incorporate  a  priori  assumptions  into 
the  reconstruction  procedure  and  examine  the  effect  these  assumptions  have  on  stability. 
Also,  given  the  continuous  dependence  result  from  Lemma  5.2  and  the  fact  that  data  is 
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invariably  noisy,  we  know  that  any  inversion  procedure  will  tend  to  give  extraneous  high 
frequency  components  in  any  estimate  of  S:  choosing  the  estimate  of  minimal  norm  should 
help  to  damp  out  spurious  components  in  the  estimate.  In  this  sense  the  procedure  may  be 
viewed  as  a  form  of  regularization. 

It  is  an  easy  application  of  Lagrange  multipliers  to  verify  that  the  unique  element  of 
with  minimum  norm  which  satisfies  the  constraints  (6.17)  must  be  of  the  form 


(6.18) 


for  some  C  C.  The  constants  A^;  can  be  determined  by  substituting  (6.18)  into 

equations  (6.17)  and  solving  the  resulting  nxn  system.  The  system  is  of  the  form  MX  =  d 
where  M  =  [m,j]  is  an  n  by  n  matrix,  A  is  the  n  vector  (Ai, . . . ,  A^)^  and  c?  is  an  n  vector 
(d(Gi), . . .  ,d{o.n)Y ■  The  entries  of  M  are  given  by 


/oo  roo 

Cj(.Ti)Cj(xi)  dxi  =  /  c(xi  —  ai)c{xi  —  Oj)  dxi. 
-OO  J  —  oo 


(6.19) 


The  matrix  M  is  clearly  Hermitian  and  in  fact  is  always  invertible.  To  see  this,  suppose  we 
can  find  some  vector  A  7^  0  with  MX  =  0.  Then  X^MX  =  0  and  we  conclude  that 

n  2 

/OO  _  TOO  ^ 

c{xi  -  o,)Aic(.ri  -  aj)Xj  dxi  =  Y 

-00  ij  J-00 

This  implies 

n 

Y  AjC(.Ti  -  aj)  =  0. 

j=i 

Fourier  transform  both  sides  and  use  the  basic  properties  of  the  Fourier  transform  to  obtain 


f(!/)  E  A, =  0. 

i=i 


(6.20) 


The  functions  fi{y)  =  are  linearly  independent  for  distinct  Oj,  and  analytic,  so  that 
f{y)  =  has  isolated  zeroes.  Based  on  equation  (6.20)  we  conclude  that  c(y)  =  0 

in  L^(IR),  contradicting  Lemma  4.2.  Therefore  M  must  be  invertible.  This  inversion  proce¬ 
dure  thus  always  produces  a  unique  estimate  of  S  if  the  measurement  locations  are  distinct. 

We  can  also  “solve”  the  inverse  problem  by  choosing  the  unique  function  S  which  satisfies 
equations  (6.17)  and  has  minimal  norm  in  a  weighted  space  T|(IR,)  with  norm  defined  by 
the  inner  product 

/•oo  _  1 

<f,g>s= 
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where  is  some  real-valued  non-negative  function  on  IR.  In  this  case,  we  have 


/oo  roo  1 

S{xi)ci{xi)dxi  =  /_^‘5(-^'i)Q(-Ti)<5(-Ti)^^d.ri 

where  we  must  assume  that  5  =  0  wherever  S  —  0.  Thus  the  integral  is  understood  to  be 
taken  only  over  that  set  where  5  is  non-zero.  Equations  (6.17)  now  take  the  form 

<S,c,5>s=d,  '  (6.21) 

and  the  minimal  norm  solution  is  of  the  form 

5(.Ti)  =  ,5(.Ti)EA,q(.ri).  (6.22) 

i=i 

The  idea  is  to  choose  (5(a'i)  to  have  the  same  general  form  as  5(.Ti),  and  so  incorporate  a 
priori  information  into  the  reconstruction  based  on  (6.22)  by  forcing  it  to  have  the  same 
general  form.  For  example,  if  we  know  that  5  is  supported  in  the  interval  [—b,  b]  we  can 
choose  S{x)  =  1  on  [— ?>,  ft]  and  (5(.t)  =  0  elsewhere.  The  optimal  estimate  of  5  becomes 

5(.t)  =  A,c,(.r) 

j=i 

where  xi-b.b]  is  the  characteristic  function  of  the  interval  [—ft,  ft]  and  where  the  A,  are  found 
by  solving 

for  j  =  1  to  n. 

Whether  we  use  a  priori  information  or  choo.se  a  uniform  weighting  on  L^,  the  recon¬ 
struction  procedure  is  as  follows:  Compute  the  matrix  M  defined  by  equation  (6.19)  and 
“measure"  the  data  d,-  at  the  corresponding  points  xi  —  a,  on  the  top  surface.  Solve  the 
system  MX  =  d  to  obtain  Ai. . . . ,  A„  and  then  compute  an  estimate  of  5(:ri)  from  equation 
(6.18)  or  (6.22).  The  stability  of  the  finite  data  inversion  is  thus  determined  by  the  nature 
of  the  matrix  M,  and  specifically,  of  its  inverse.  We  can  quantify  the  stability  of  the  finite 
data  inverse  problem  by  studying  the  conditioning  of  the  inverse  of  M  in  various  situations. 
This  is  done  in  the  following  section  by  studying  the  singular  values  of  the  matrix  M. 
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7  Numerical  Experiments 

We  will  now  examine  the  finite  data  version  of  the  inverse  problem  using  the  previously 
described  inversion  procedure.  In  this  section  we  apply  the  procedure  to  simulated  data 
sets,  both  with  and  without  noise.  Our  main  focus  is  to  examine  the  stability  and  resolution 
of  back  surface  estimates  with  respect  to  various  experimental  parameters,  specifically,  the 
frequency  of  the  input  heat  flux  and  the  distribution  of  the  measurement  locations  along 
the  top  surface  of  the  sample.  We  also  demonstrate  how  a  priori  assumptions  about  the 
nature  of  the  corrosion  can  be  incorporated  into  the  inversion  procedure,  and  the  effects 
such  assumptions  have  on  stability  and  resolution. 

There  are  a  few  points  worth  mentioning  before  we  present  the  numerical  results.  The 
estimates  of  S  constructed  using  equation  (6.18)  lie  in  H°^{WC)  and  so  the  graph  X2  =  S{xi) 
would  make  sense  as  a  curve  in  except  that  the  estimate  will  usually  be  complex- valued. 
This  is  not  surprising,  for  a  complex-valued  S  makes  perfect  sense  in  the  linearized  version 
of  the  prolDlem  on  which  the  inversion  procedure  is  based.  Estimates  of  S,  especially  in  the 
presence  of  noise  or  the  linearization  error,  will  almost  certainly  have  non-zero  imaginary 
part.  However,  for  those  S  of  “small”  norm  the  linearized  problem  accurately  reflects  the 
full  non-linear  (with  respect  to  S)  direct  problem  and  so  the  estimate  of  S  should  be  “mostly 
real” ,  that  is,  it  should  have  a  relatively  small  imaginary  component.  This  is  indeed  the  case. 
A  physically  meaningful  estimate  of  the  true  back  surface  X2  =  ^(-Ti)  can  be  provided  by 
either  dropping  the  imaginary  component  or  taking  the  modulus  of  the  estimate.  We  choose 
the  latter. 

In  the  examples  that  follow  we  generate  simulated  test  data  using  the  full  direct  problem 
(2.2)  with  heating  g{x)  =  1.  The  direct  problem  is  solved  by  converting  it  into  a  boundary 
integral  equation  which  is  then  solved  numerically.  The  boundary  integral  formulation  leads 
to  a  second  kind  Fredholm  equation  (explained  below).  The  fact  that  g  is  not  compactly 
supported  nor  even  presents  a  minor  problem.  This  can  be  fixed  by  simply  subtracting 
off  the  function  uq  satisfying  the  direct  problem  with  5  =  0.  Recall  that  uq  can  be  found 
explicitly  by  solving  equation  (3.5).  The  function  v  =  ii  —  hq  then  satisfies  the  boundary 
value  problem 
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Av  —  iwv  =  0  in  Q, 
dv 

—  —  0  on  .T2  =  1, 
av 

dv  diiQ 


—  =  on  .T2  =  5(.ri). 

ov  av 

If  S'(.ri)  is  sufficiently  rapidly  decreasing  then  the  Neumann  data  in  the  above  boundary 
value  problem  is  and  the  integral  equation  formulation  becomes 

-;^u(.r)  +  f  ^(:r,  y)v{y)  day  ^  [  G(:r,  y)^  day 

2  JdQo  avy  Jano  ov 

where  day  is  surface  measure  on  dflo  and  G{x,y)  =  r(|.r  —  ?/|).  The  integral  equation  is 
solved  using  Nystrom’s  method  (see  [2])  with  appropriate  quadrature  rules  on  the  interval 
(—00,00)  that  arises  when  the  boundary  of  Qq  is  parameterized.  The  boundary  integral 
approach  is  exceptionally  fast  and  accurate  and  provides  estimates  of  the  solution  only  on 
the  boundary  of  flo.  the  only  place  that  we  need  the  solution. 

To  illustrate  the  general  procedure  and  to  show  that  the  inversion  algorithm  provides 
reasonable  estimates,  we  begin  with  a  simple  example.  We  apply  the  inversion  procedure  to 
data  generated  using  the  back  surface 


Six)  = 


,-(t+2)V2 


We  use  a  heating  frequency  of  w  =  1.  Asa  first  step  the  function  c(.t)  defined  in  Lemma  4.1 
is  computed  by  solving  the  system  (4.10)  with  a  =  0  numerically  and  then  computing 
d)o  =  — 2r  -f  uo,  again  via  Nystrom’s  method.  These  quantities  are  precomputed  and  stored 
rather  than  re-computed  every  time  they  are  needed,  for  the  function  c(;r)  depends  only  on 
flo.  not  S.  The  same  is  true  of  the  matrix  M.  The  entries  of  the  matrix  are  computed 
by  using  equation  (6.19)  and  an  appropriate  quadrature  rule  on  (—00,00).  In  the  example 
below  the  data  is  computed  using  the  full  direct  problem.  The  temperature  data  vector  d 
is  computed  at  21  equally  spaced  points  on  the  top  surface,  .Ti  =  0,  where  o,  =  — 5  -f  |  for 
i  =  0  to  20.  We  then  invert  the  21  x  21  s.ysteni  MX  =  d  to  find  A  and  return  an  estimate  of 
S  via  equation  (6.18).  The  estimate  of  S  is  computed  at  a  suitable  number  of  points  on  the 
range  of  interest,  in  this  case  from  —5  to  5.  The  reconstruction  is  shown  in  Figure  2.  The 
dotted  line  is  the  actual  function  5(.t)  and  the  solid  line  is  the  reconstructed  version. 
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Stability 

Of  particular  interest  is  the  sensitivity  of  the  inversion  procedure  with  respect  to  various 
experimental  parameters,  e.g.,  heating  frequency  and  measurement  locations.  The  first  task 
is  to  quantify  the  stability  or  conditioning  of  the  finite  data  inverse  problem.  One  sensible 
way  to  do  this  is  to  perform  a  singular  value  decomposition  on  the  matrix  M  defined  by 
equation  (6.19)  and  examine  the  magnitude  of  the  singular  values.  When  the  singular  values 
are  small  the  inversion  of  MX  =  d  magnifies  small  perturbations  in  d.  Put  another  way, 
small  singular  values  mean  that  relatively  large  changes  in  S  (and  so  in  A)  produce  relatively 
small  changes  in  the  data,  so  that  perturbations  in  the  back  surface  are  “hard  to  see.”  Our 
goal  in  choosing  experimental  parameters  is  therefore  to  make  the  singular  values  of  M  as 
large  as  possible,  within  certain  limits.  We  should  remark  that  one  can  define  the  condition 
number  of  M  as  the  ratio  of  the  largest  to  smallest  singular  values  and  then  attempt  to 
quantify  stability  using  this  single  number.  This  is  not  always  good  approach  in  the  present 
setting,  as  later  examples  will  show. 

Let  us  begin  by  examining  how  the  stability  of  the  inversion  procedure  depends  on  the 
locations  of  the  temperature  measurements  on  the  top  surface.  In  the  following  examples 
we  fix  the  heating  frequency  at  cj  =  1  and  take  measurements  of  the  resulting  temperature 
at  21  equally  spaced  locations  on  the  interval  [—a,  a]  for  several  values  of  a.  The  resulting 
measurement  locations  are  therefore  of  the  form  a,-  =  — o  +  ^  for  ?'  =  0, . . . ,  20.  In  each  case 
the  matrix  M  is  computed  and  a  singular  value  decomposition  is  performed.  Let  the 
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—  —  —  a  =  10.0 

-  a  =  5.0 

- a  =  3.0 

- a  =  2.0 

- a  =  1.0 


Figure  3:  logjo|Ot  |  versus  i  for  various  values  of  a. 

singular  values  of  M  be  denoted  by  q,-,  i  =  1  to  21,  arranged  in  descending  order.  In  Figure  3 
we  plot  the  quantity  logjo  versus  i  for  the  cases  a  =  1, 2, 3, 5, 10. 

It  is  apparent  that  as  the  measurement  locations  become  more  spread  out  (as  a  gets 
larger)  the  singular  values  decay  more  slowly  and  hence  the  inversion  procedure  becomes 
more  stable.  In  light  of  Theorem  5.2  this  is  not  surprising.  When  the  measurement  locations 
are  close  together  we  are  able  to  resolve  higher  spatial  frequencies  in  the  data  and  so  we 
are  able  to  estimate  higher  frequencies  in  the  Fourier  decomposition  of  S.  But  according  to 
Theorem  5.2,  these  are  exactly  the  portions  of  S  that  are  difficult  to  reconstruct — they  are 
heavily  damped  out  in  the  data.  The  finite  data  version  of  the  problem  reflects  this,  with  a 
full  6  orders  of  magnitude  variation  for  the  smallest  singular  values  between  the  cases  a  =  1 
and  a  =  10. 

Another  way  to  look  at  the  stability  of  the  various  experimental  configurations  is  to  sup¬ 
pose  that  we  have  an  “error  magnification  tolerance’'  E,  and  that  in  the  inversion  procedure 
we  disregard  all  singular  vectors  whose  singular  values  are  less  than  This  idea  has  been 
used  in  studying  the  stability  for  the  impedance  imaging  problem  [8].  The  inversion  proce¬ 
dure  is  then  stabilized  at  the  expense  of  rendering  those  components  of  S  lying  in  the  span 
of  the  corresponding  functions  invisible.  Figure  4  shows  the  number  of  singular  values  of  M 
which  satisfy  ^  versus  logio(F^)  for  E  from  1  to  10“®.  As  in  the  previous  examples,  the 
matrix  M  is  21  x  21  and  we  use  measurement  locations  on  the  top  surface  o,-  =  —  a  -f 
?'  =  0, . . . ,  20  for  a  =  1, 2, 3, 5, 10.  The  heating  frequency  is  a’  =  1. 
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Figure  4-'  Number  of  singular  values  with  aj  >  ^  versus  logio(£')  for  various  values  of  a. 

Figure  4  also  makes  clear  that  as  the  measurement  locations  become  spread  out  more 
singular  values  satisfy  a,  >  ^.  The  inversion  procedure  then  admits  more  basis  functions, 
presumably  improving  the  fidelity  of  the  reconstruction.  In  the  two  cases  below  we  perform 
the  actual  reconstruction  with  E  =  100  (so  only  singular  values  greater  than  0.01  are  ad¬ 
missible)  and  add  a  small  amount  of  random  noise  to  the  data  (equal  to  10  percent  of  the 
maximum  signal  strength).  We  then  perform  a  reconstruction  which  omits  all  basis  vectors 
whose  corresponding  singular  values  are  less  than  Figure  5  illustrates  the  case  in  which 
the  measurements  locations  are  equally  spaced  from  —5  to  5;  there  are  9  admissible  singular 
values. 


Figure  5:  Reconstruction  of  S{x)  for  21  measurements  on  [—5,5],  tolerance  E  =  10^. 

In  Figure  6  we  take  the  21  measurements  on  the  smaller  interval  [—1, 1],  which  yields  only  3 
admissible  singular  values. 
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Figure  6:  Reconstruction  of  S{a-)  for  21  measurements  on  [—1, 1],  tolerance  E  =  10^. 

The  reconstruction  in  Figure  6  is  noticeably  inferior  to  that  of  Figure  5,  but  we  have  only  3 
admissible  basis  functions  with  which  to  construct  S{x).  Increasing  the  value  of  E  to  admit 
more  basis  functions  is  not  successful.  Figure  7  illustrates  what  happens  if  we  take  jF  =  10^* 
with  measurements  on  [—1,1].  Now  5  singular  values  are  admissible,  but  the  reconstruction 
is  overwhelmed  by  noise. 


Figure  1:  Reconstruction  of  S(x)  for  21  measurements  on  [—1,1],  tolerance  E  =  10^. 

The  moral  seems  clear;  for  maximum  stability  with  a  fixed  number  of  measurement 
locations,  we  should  spread  the  measurements  over  as  large  a  region  as  possible.  There  are 
limits  to  this  approach,  however.  If  we  spread  out  the  measurements  we  do  gain  stability, 
but  we  will  no  longer  be  able  to  estimate  high  frequencies  in  the  Fourier  decomposition  of 
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S.  This  is  illustrated  by  Figure  8,  where  we  take  21  noise-free  measurements  on  the  interval 
[-10, 10]  and  estimate  5  with  error  tolerance  E  =  10^.  In  this  case  all  of  the  singular  values 
are  admissible. 


Figure  8:  Reconstruction  of  S{x)  for  21  measurements  on  [-10,10],  tolerance  E  =  10^. 

Despite  the  fact  that  the  inversion  is  quite  stable,  our  inability  to  resolve  high  frequencies 
results  in  a  loss  of  resolution  of  small-scale  detail  in  the  reconstruction.  With  regard  to  the 
distribution  of  the  measurement  locations,  the  reconstruction  process  involves  a  compromise 
between  stability  and  resolution  of  small-scale  features. 

In  the  next  series  of  examples  we  examine  the  dependence  of  the  stability  on  oj,  the 
frequency  of  the  applied  heat  flux.  We  consider  the  cases  oo  =  0.01,0.1,1.0,10.0,100.0. 
In  each  case  we  take  21  equally  spaced  temperature  measurements  on  the  interval  [—5,5]. 
Figure  9  shows  log^g  |ai|  versus  i  for  each  case. 


w  =  0.01 
w  =  0.1 
w=  1.0 
w=10.0 
w=  100.0 


Figure  9:  log^g  |o?|  versus  i  for  various  values  of  u. 

The  figure  illustrates  that  higher  frequencies  give  rise  to  much  smaller  singular  values. 
As  before,  it  is  instructive  to  consider  the  case  in  which  we  have  an  error  tolerance  E  and  in 
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the  inversion  process  we  omit  those  singular  vectors  whose  singular  values  are  less  than 
Figure  10  shows  the  number  of  singular  values  which  exceed  ^  for  E  from  10“*^  to  10^’’,  for 
each  of  the  frequencies. 


w  =  0.01 
w  =  0.1 
w  =  1.0 
w  =  10.0 
w  =  100.0 


Figure  10:  Number  of  singular  values  with  o,  >  ^  versus  logio(£')  for  various  values  ofco. 


Figure  10  clearly  illustrates  the  situation.  For  E  =  10~^  as  before,  a;  =  10  and  u  =  100 
have  no  admissible  singular  values  at  all.  The  reconstruction  (if  carried  out)  is  identically 
zero.  The  smaller  singular  values  at  higher  frequencies  are  due  to  the  fact  that  at  higher 
frequencies  the  periodic  heating  penetrates  very  little  into  the  sample  and  becomes  more  of 
a  “skin  effect."  As  a  result  very  little  energy  reaches  the  back  surface  and  even  less  returns 
to  be  measured  on  the  top  surface;  at  high  frequencies  the  map  taking  S  into  d  is  essentially 
multiplication  by  zero.  It  is  interesting  to  note  that  while  higher  frequencies  produce  smaller 
singular  values,  the  condition  number  of  M  for  oj  =  100  is  only  10.3,  while  the  condition 
number  for  uj  =  0.01  is  706.1.  Clearly,  though,  it's  not  enough  to  make  the  condition  number 
small.  The  singular  values  themselves  must  be  large  enough  for  the  inversion  procedure  to 
be  stable  in  the  presence  of  a  fixed  noise  level. 

While  lower  frequencies  make  the  inversion  process  more  stable,  there  are  limits  to  how 
small  we  can  make  u  and  maintain  resolution.  Figure  11  illustrates  a  reconstruction  based 
on  u!  =  0.1.  The  parameters  are  otherwise  identical  to  those  that  were  used  to  produce 
Figure  5.  All  of  the  singular  values  are  admissible.  In  fact,  the  smallest  singular  value  is 
0.063.  As  with  the  case  in  which  the  measurement  locations  were  spread  out  over  [—10, 10], 
we  lose  resolution  at  low  temporal  frequencies  for  the  input  flux. 
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Figure  11:  Reconstruction  of  S{x)  for  21  measurements  on  [—5,5],  a;  =  0.1,  tolerance 

E  =  102. 

Incorporating  A  Priori  Information 

The  preceding  examples  illustrate  that  the  inversion  procedure  involves  a  compromise 
between  stability  and  resolution.  If  the  data  points  are  too  closely  spaced,  the  inversion 
procedure  is  unstable.  If  the  data  points  are  too  spread  out,  the  inversion  procedure  becomes 
stable,  but  resolution  is  lost;  measurements  taken  far  from  the  support  of  the  defect  contain 
little  information,  because  the  heat  diffuses  veiy  rapidly.  Variations  in  the  input  heating 
frequency  give  rise  to  a  similar  phenomena.  How  shall  we  find  the  ‘"best”  experimental 
parameters?  One  useful  possibility  is  to  incorporate  a  priori  information  or  assumptions  into 
the  inversion  procedure.  We  will  illustrate  the  idea  by  examining  the  problem  under  the 
assumption  that  the  defect  or  function  S  is  supported  in  a  known  interval. 

In  the  following  examples  we  assume  that  the  defect  being  imaged  is  supported  in  the 
interval  [—2,2].  The  only  modification  to  the  inversion  procedure  is  that  the  matrix  M  is 
computed  in  accordance  with  equation  (6.21)  and  the  function  S  is  estimated  using  equation 
(6.22).  We  will  study  the  stability  of  the  inversion  procedure  with  respect  to  the  distribution 
of  the  measurement  locations  on  the  top  surface. 

As  in  the  previous  cases,  we  choose  measurement  locations  at  Xi  =  a,  on  the  sample  top 
surface,  where  Oj  =  — a  +  for  ?  =  0  to  20.  The  heating  frequency  in  all  cases  that  follow 
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a  =  0.5 
a=  1.0 
a  =  2.0 
a  =  5.0 
a  =  10.0 


Figure  12:  Singular  values  o,  versus  i  for  various  values  of  a. 

is  a;  =  1.  Let  us  begin  by  examining  the  singular  values  of  the  inversion  matrix  AI  for  a  few 
choices  of  a.  In  Figure  12  we  plot  the  quantity  logjo  |a,  |  versus  i  for  a  =  0.5, 1.0, 2.0, 5.0, 10.0. 

The  figure  shows  that  the  best  conditioning  for  the  inverse  problem  occurs  at  o  =  2, 
when  the  measurement  locations  are  distributed  approximately  in  the  same  interval  in  which 
the  defect  is  assumed  to  be  supported.  As  before,  closely  spaced  locations  give  rise  to 
an  ill-conditioned  problem.  However  unlike  the  previous  cases  widely  spaced  nodes  also 
result  in  poor  conditioning.  When  AI  is  computed  using  equation  (6.21)  those  rows  of  M 
corresponding  to  measurement  locations  far  from  the  support  of  S  are  very  nearly  set  to  zero 
since  the  function  c(;r  —  o,  )  is  rapidly  decreasing  away  from  Oj. 

If  an  error  magnification  tolerance  E  is  specified,  we  can  plot  the  number  of  allowable 
singular  values  o,  >  ^  versus  logio(F)  for  the  different  node  spacings: 

-  a  =  0.5 

—  —  —  a  =  1 .0 

-  a  =  2.0 

- .  a  =  5.0 

- a  =10.0 


Figure  13:  Number  of  singular  values  with  o,  >  ^  versus  logjolIF)  for  various  values  of  a. 

As  expected,  a  =  2.0  allows  more  singular  values  for  a  fixed  value  of  E  than  any  other 
choice  for  measurement  spacing.  It  is  useful  to  look  at  a  few  reconstructions  based  on  this 
strategy.  In  the  two  cases  below  we  take  E  —  300  (so  only  singular  values  greater  than  ^ 
are  admissible)  and  add  a  small  amount  of  random  noise  to  the  data  (equal  to  1  percent  of 
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the  maximum  signal  strength).  We  then  perform  a  reconstruction  which  omits  all  singular 
values  less  than  The  function  defining  the  back  surface  is  S{x)  = 

Figure  14  illustrates  the  first  case  using  a  =  2,  the  best  choice  according  to  Figure  13.  In 
this  case  7  singular  values  are  admissible. 


Figure  14:  Reconstruction  of  S{x)  for  21  measurements  on  [-2,2],  tolerance  E  =  300. 
For  a  -  10  we  have  4  admissible  singular  values  and  the  reconstruction  shown  in  Figure  15. 


Figure  15:  Reconstruction  of  S{x)  for  21  measurements  on  [—10,10],  tolerance  E  =  300. 

The  case  a  =  0.5  also  yields  4  admissible  singular  values  and  the  reconstruction  shown  in 
Figure  16. 
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Figure  16:  Reconstruction  of  S{x)  for  21  measurements  on  [—0.5, 0.5],  tolerance  E  =  300. 

The  actual  reconstructions  confirm  that  a  =  2  yields  the  most  desirable  results.  Choosing 
a  significantly  smaller  or  larger  than  the  support  of  S  results  in  decreased  stability  and/or 
accuracy  for  the  reconstruction. 

Of  course,  the  assumption  that  S  is  supported  in  a  given  interval  should  be  detrimental 
to  the  reconstruction  if  that  assumption  turns  out  to  be  false.  In  the  following  case  we  let 

and  perform  the  reconstruction  under  the  assumption  that 
S  is  supported  in  the  interval  [—2,2].  We  take  measurements  at  21  equally  spaced  location 
between  —2  and  2,  the  best  case  from  above,  and  use  an  error  tolerance  E  =  300.  The  result 


Figure  17:  Reconstruction  of  S{x)  for  21  measurements  on  [—2,2],  tolerance  E  =  300. 
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The  incorrect  assumption  obviously  introduces  errors  into  the  reconstruction,  although  that 
portion  of  S  which  is  non-zero  in  the  interval  [—2, 2]  is  still  recovered  with  reasonable  accu¬ 
racy. 

Concluding  Remarks. 

In  this  paper  we  have  investigated  the  inverse  problem  of  recovering. an  unknown  bound¬ 
ary  portion  of  some  object  by  applying  a  heat  flux  to  an  accessible  portion  of  the  boundary 
and  measuring  the  resulting  temperature  response.  We  have  considered  a  linearized  version 
of  the  problem  and  found  that  the  continuous  version  of  the  inverse  problem,  in  which  one 
has  data  at  every  point  on  the  accessible  portion  of  the  surface,  is  extremely  ill-posed.  In¬ 
deed,  the  linearized  version  requires  one  to  solve  a  first  kind  convolution  integral  equation 
for  the  unknown  surface.  The  convolution  kernel  has  a  Fourier  transform  which  dies  rapidly 
at  infinity,  and  so  the  inversion  is  extremely  sensitive  to  the  data  at  high  spatial  frequencies. 
We  performed  a  variety  of  numerical  studies  which  show  that  the  ill-posedness  is  directly 
reflected  in  the  finite  data  version  of  the  problem,  by  the  rapid  decay  of  the  singular  values 
of  the  matrix  which  governs  the  inversion  process.  This  ill-posedness  depends  on  a  num¬ 
ber  of  factors;  in  particular,  the  locations  of  the  measurements  have  a  large  effect  on  the 
conditioning  of  the  inverse  problem,  and  these  effects  mirror  the  behavior  of  the  continuous 
version.  We  have  also  considered  the  effect  of  including  a  priori  assumptions  in  the  finite 
data  inversion  procedure,  by  weighting  appropriate  Hilbert  spaces  in  which  the  solution  S 
resides.  The  inclusion  of  this  information  can  help  in  determining  the  optimal  locations  for 
measurements  on  the  top  surface. 

There  are  a  number  of  interesting  directions  we  could  take  from  here.  In  our  studies 
we  used  only  the  input  flux  whose  magnitude  is  identically  one  on  the  top  surface.  Similar 
results  can  l^e  obtained  for  more  general  fluxes,  and  this  would  allow  one  to  study  the  effect 
that  the  input  heat  flux  has  on  sensitivity  and  resolution.  The  fully  time-dependent  case 
would  also  be  of  interest.  The  procedure  presented  in  this  paper  would  also  work  for  a  full 
three-  dimensional  problem,  although  qualitatively  the  results  should  be  the  same — the  high 
spatial  frequencies  in  the  back  surface  should  be  difficult  to  see. 

As  mentioned  earlier,  the  inversion  process  which  chooses  that  function  with  minimal 
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norm  which  is  consistent  with  the  measured  data  seems  to  act  like  a  form  of  regularization  for 
the  inverse  problem.  It  would  be  interesting  to  examine  in  what  sense  this  is  true,  and  how 
it  relates  to  more  traditional  forms  of  regidarization.  It  is  also  possible  (and  not  difficult)  to 
carry  out  the  same  minimization  process  in  higher  Sobolev  spaces,  e.g.,  ,  and  thus  put  a 

higher  “penalty”  on  functions  with  oscillations.  This  too  would  make  an  interesting  study. 
We  would  also  like  to  examine  conditions  under  which  our  inversion  procedure  is  guaranteed 
to  converge  to  the  solution  of  the  linearized  inverse  problem. 
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